#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Feb  9 14:46:52 2021

@author: liqingsimac
"""

'''
# 4.1.1.
import numpy as np
xc,xd = np.linspace(0,1,11,retstep=True)
xo = np.linspace(0,1,10,endpoint=False)
xlog = np.logspace(0,100,2)
x = np.arange(1,10)
z = np.zeros(10, dtype=float)
y = np.ones(10, dtype=float)
la = np.array([1,2,3,0])
xe = np.empty_like(x)
'''

'''
# 4.1.2.
import numpy as np
x = np.linspace(0,1,11)
y = np.zeros_like(x)

xL=np.linspace(0, 0.5, 5, endpoint=False)
xM=np.linspace(0.5, 0.6, 10, endpoint=False)
xR=np.linspace(0.6, 1.0, 5)
xs=np.hstack( (xL,xM,xR) )
'''

'''
# 4.1.3. 
import numpy as np
a = np.linspace(0,1,5)
c = np.linspace(1,3,5)
a+c, a*c, a/c

import numpy as np
a = np.linspace(0,1,5)
2*a, a*2, a/5, a**3, a+2

a=np.ones(4,dtype=int)
b=np.linspace(0,1,4)
a+=b
'''


'''
# 4.1.3. 
import numpy as np
f=np.array([1,1,0,0,1,2])

f_av=f.copy()

for i in range(1,len(f)-1):
    f_av[i]=(f[i-1]+f[i]+f[i+1])/3.0

f_ave=f.copy()
f_ave[1:-1]=(f[:-2]+f[1:-1]+f[2:])/3.0
'''


'''
# 4.1.4. 
import numpy as np
def h(x):
#     return 1 if 0<=x<=1 else 0 
    if x<0.0:
        return 0.0
    elif x<=1.0:
        return 1.0
    else:
        return 0.0
    
v=np.linspace(-2,2,41)
hv=h(v)
'''


'''
# 4.1.3. 
a=[1,2,3,4,5]
b=a[1:-1]
b
b[1]=100
b
a
x=np.array([1,2,3,4,5])
x
y=x[1:-1]
y
y[1]=100
y
x
'''


'''
# 4.1.5. 
import numpy as np
x=np.linspace(-2,2,9)
y=x<0

z=x.copy()
z[y]=-z[y]
z

w=x.copy()
w[w<0]=-w[w<0]
w

h=np.ones_like(x)
h[x<0]=0
h[x>1]=0

'''


'''
# 4.1.5. 
import numpy as np
import matplotlib.pyplot as plt
x=np.linspace(-1,3,9)
choices=[ x>=2, x>=1, x>=0, x<0 ]
outcomes=[ 4.0, x**2, x**3, -x ]
k=np.select(choices,outcomes)
plt.plot(x,k,'-o')
plt.savefig('pic/fig-4-1-5.png')


import numpy as np
import matplotlib.pyplot as plt
def m1(x):
    return np.exp(2*x)
def m2(x):
    return 1.0
def m3(x):
    return np.exp(1-x)

x = np.linspace(-10,10,21)
conditions = [ x>=0, x>=1, x<0 ]
functions = [ m2, m3, m1 ]
m = np.piecewise(x, conditions, functions)
plt.plot(x,m,'-o')
plt.savefig('pic/fig-4-1-5-piecewise.png')
'''


'''
# 4.2. 
import numpy as np
x=np.linspace(0,1.0,11)
x.ndim
x.shape
x.dtype

x=np.array([[0,1,2,3],[10,11,12,13],[20,21,22,23]])
x.ndim
x.shape
x.dtype
'''


'''
# 4.2.1. 
x=np.array([[0,1,2,3],[10,11,12,13],[20,21,22,23]])
r=np.array([2,3,4,5])
c=np.array([[5],[6],[7]])

print(2*x)
print(r*x)
print(x*r)
print(c*x)
'''


'''
# 4.2.2. 
xv = np.linspace(-1,1,5)
yv = np.linspace(-1,1,3)
[xa,ya] = np.meshgrid(xv,yv)
print('xv=\n',xv)
print('yv=\n',yv)
print('xa=\n',xa)
print('ya=\n',ya)
print('xa*ya=\n',xa*ya)
'''


'''
# 4.2.2. 
x1=np.mgrid[-1:1:9j]
print(x1)
x2=np.mgrid[-1:1:0.25]
print(x2)

[xm,ym]=np.mgrid[-1:1:5j, 0:1:3j]
print('xm=\n',xm)
print('ym=\n',ym)
print('xm*ym=\n',xm*ym)

[xo,yo] = np.ogrid[-1:1:5j, 0:1:3j]
xo
yo

x=np.zeros((4,3), dtype=float)
'''

'''
# 4.2.3. 
x=range(6)
a=np.reshape(x,(2,3))
print('a=',a)
print('the type of a is',type(a))

v=np.linspace(0,1,5)
print('v=',v)
vg=np.reshape(v,(5,1))
print('the shape of vg is',vg.shape)
print('vg=',vg)
'''


'''
# 4.2.4. 
u=np.linspace(100,300,3)
vg=np.reshape(np.linspace(0,28,15),(5,3))
print('u=',u)
print('vg=',vg)
print('u+vg=',u+vg)
print('a slicing of vg is:',vg[1:-1,1:])
'''


'''
# 4.3. 
[xo,yo,zo]=np.ogrid[1:2:2j,30:50:3j,600:900:4j]
print('xo=',xo)
print('yo=',yo)
print('zo=',zo)
print('xo+yo+zo=',xo+yo+zo)
'''


'''
# 4.4.1.
import numpy as np
quarter=np.array([1,2,3,4],dtype=int)
results=np.array([37.4, 47.3, 73.4, 99])
outfile=open('q4.txt','w')
outfile.write('The results for the first four quarters\n\n')
for q,r in zip(quarter,results):
    outfile.write('For quarter %d the result is %5.1f\n'%(q,r))
outfile.close()
'''


'''
# 4.4.1. 
infile=open('q4.txt','r')
lquarter=[]; lresult=[]
temp00=infile.readline(); temp01=infile.readline()
for line in infile:
    words=line.split()
    lquarter.append(int(words[2]))
    lresult.append(float(words[6]))
infile.close()

import numpy as np
aquarter=np.array(lquarter,dtype=int)
aresult=np.array(lresult)
print('quarters=',aquarter)
print('results=',aresult)
'''


'''
# 4.4.2. 
import numpy as np
leng=21
x=np.linspace(0,2*np.pi,leng)
c=np.cos(x)
s=np.sin(x)
t=np.tan(x)
arr=np.empty((4,leng),dtype=float)
arr[0,: ]=x
arr[1,: ]=c
arr[2,: ]=s
arr[3,: ]=t
np.savetxt('x.txt',x)
np.savetxt('xcst.txt',(x,c,s,t))
np.savetxt('xarr.txt',arr)


xc=np.loadtxt('x.txt')
xc,cc,sc,tc=np.loadtxt('xcst.txt')
arrc=np.loadtxt('xarr.txt')
'''


'''
# 4.4.3. 
np.save('array.npy',arr)
arrd=np.load('array.npy')

np.savez('test.npz',x=x,c=c,s=s,t=t)
temp=np.load('test.npz')
print(temp.files)
xc=temp['x']; print('xc=',xc)
cc=temp['c']
sc=temp['s']
tc=temp['t']
'''


'''
# 4.6. 
import numpy as np
x=np.array([[5,4,1],[7,3,2]])
np.max(x)
np.max(x,axis=0)
np.max(x,axis=1)
'''

'''
# 4.7.2. 
import numpy as np
roots=[0,1,1,2]
coeffs=np.poly(roots)
print('coeffs=',coeffs)
'''


'''
# 4.7.2. 
import numpy as np
x=np.linspace(0,2*np.pi,21)
y=np.sin(x)
c=np.polyfit(x,y,3)
print('c=',c)

import matplotlib.pyplot as plt
y1=np.polyval(c,x)
plt.plot(x,y,'bo-',x,y1,'r--')
plt.legend(['y=sin(x)','polynomial fit'])
plt.savefig('pic/fig-4-7-2.png')
'''


'''
# 4.8. 
import numpy as np
A=np.matrix([[1,2],[3,4]])
b=np.matrix([[5],[6]])
print(A*b)

A=np.array([[1,2],[3,4]])
b=np.array([[5],[6]])
print(np.dot(A,b))

E=np.identity(3,dtype=float)
C=2*np.eye(3,4,-1)+3*np.eye(3,4,0)+4*np.eye(3,4,1)

v1=np.array([1,2,3])
v2=np.array([4,5,6])
rows=np.vstack((v1,v2))
print(rows)
print(rows.T)
'''


'''
# 4.8.2. 
import numpy as np
import numpy.linalg as npl

A=np.array([[4,2,0],[9,3,7],[1,2,1]])
print('A=',A)
print('det(A)=',npl.det(A))
B=npl.inv(A)
print('The inverse of A is: ',B)
print('B*A=',np.dot(B,A))


import numpy as np
import numpy.linalg as npl
A=np.array([[-2,-4,2],[-2,1,2],[4,2,5]])
evals,evecs = npl.eig(A)
eval1=evals[0]
evec1=evecs[:,0]
print(eval1)
print(evec1)
'''


'''
# 4.8.3. 
import numpy as np
import numpy.linalg as npl
A=np.array([[3,2,1],[5,5,5],[1,4,6]])
b=np.array([[5,1],[5,0],[-3,-7/2]])
x=npl.solve(A,b)
print(x)
np.dot(A,x)-b
'''


'''
# 4.9.1. 
import numpy as np
import scipy.optimize as sco

def fun(x):
    return np.cosh(x)/np.sinh(x)-x

roots=sco.fsolve(fun,1.0)
root=roots[0]
print('The root is: %15.12f' % root)
print('The value is: %e' % fun(root))


x=np.linspace(-3,3,21)
y1=x
y2=np.cosh(x)/np.sinh(x)

import matplotlib.pyplot as plt
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(x,y1,'b-',x,y2,'r-')
ax.legend(['y=x','y=coth(x)'])
ax.axis('equal')
ax.set(xlim=(-3,3),ylim=(-2,2))
ax.set(xlabel='x')
ax.set(ylabel='y')
ax.axvline(x=0)
ax.axhline(y=0)
ax.set_title('Solve coth(x)=x')
fig.savefig('pic/fig-4-9-1.png')
'''


